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We study the prospects for three-dimensional mapping of the dark matter to high redshift through 
the shearing of faint galaxies images at multiple distances by gravitational lensing. Such maps could 
provide invaluable information on the nature of the dark energy and dark matter. While in principle 
well-posed, mapping by direct inversion introduces exceedingly large, but usefully correlated noise 
into the reconstruction. By carefully propagating the noise covariance, we show that lensing contains 
substantial information, both direct and statistical, on the large-scale radial evolution of the density 
field. This information can be efficiently distilled into low-order signal-to-noise eigenmodes which 
may be used to compress the data by over an order of magnitude. Such compression will be useful 
for the statistical analysis of future large data sets. The reconstructed map also contains useful 
information on the localization of individual massive dark matter halos, and hence the dark energy 
from halo number counts, but its extraction depends strongly on prior assumptions. We outline 
a procedure for maximum entropy and point-source regularization of the maps that can identify 
alternate reconstructions. 



I. INTRODUCTION 

Embedded within the evolution of the three- 
dimensional distribution of the dark matter lies a wealth 
of information on the nature of the dark energy and dark 
matter in the universe. The growth of its clustering in 
volumes associated with redshift, along with the abun- 
dance of discrete dark matter clumps or halos which it 
controls, is one of our most direct probes of the expan- 
sion history (e.g. Q). It is certainly the one that is best 
understood from the theoretical standpoint. 

Unfortunately, most probes of structure in the low red- 
shift universe rely on luminous matter, e.g. galaxies and 
clusters of galaxies, as tracers of dark matter distribu- 
tion. Such probes are subject to significant uncertainties 
in the physical processes that govern the formation and 
evolution of the tracers. The only exception is the image 
distortion from gravitational lensing of distant objects 
by the dark matter. On large scales and for small distor- 
tions, this is known as weak gravitational lensing j^, [|. 
The distortion of faint galaxy images by the large-scale 
structure of the universe has now been detected with high 
significance by several experimental groups [Q. 

A fundamental obstacle for weak lensing studies of the 
matter distribution is that the technique is inherently 
two-dimensional. All of the matter along the line-of-sight 
to a distant source contributes to lensing and so the dis- 
tortion reflects a two-dimensional projection of the dark 
matter. Unfortunately then the evolution of structure 
is hidden in the missing radial dimension. This limita- 
tion can in principle be overcome by a tomographic re- 
construction of the three-dimensional distribution from 
sources spanning a range of distances or redshifts. 

Under fairly restrictive assumptions, this tomographic 
technique has been applied to lensing data to localize 



the halo associated with a cluster of galaxies jj| and val- 
idated by follow-up studies. The critical assumption is 
that the lensing mass be a single halo, well localized in 
redshift. Taylor Q has recently shown that these and 
other restrictions are unnecessary in principle. In the 
absence of noise, tomographic mapping of the dark mat- 
ter is a well-posed problem. In this paper, we study the 
feasibility of reconstructing three-dimensional dark mat- 
ter maps in the presence of noise. We will show that 
a careful accounting of the noise and in particular its 
covariance across the map is essential for extracting in- 
formation from the map. 

Tomography also presents a severe data analysis chal- 
lenge, similar to but potentially far worse than that fac- 
ing cosmic microwave background (CMB) experiments. 
A full weak lensing data set will have a two-component, 
two-dimensional megapixel map for each of ten or more 
source redshift slices. We also study how techniques de- 
veloped for the CMB and galaxy redshift surveys may 
be applied to compress these data to a more manageable 
size. 

The outline of the paper is as follows. We begin in 
§ II with a brief review of mapmaking techniques and 
apply them to the two-dimensional lensing observables 
keeping careful track of the propagation of measurement 
errors. We use these techniques to reconstruct the three- 
dimensional distribution of the dark matter in 8 [II. Al- 



though the reconstruction is extremely noisy for any re- 
alistic situation, the noise has very particular properties 
that are absent in the signal. This fact is used to regu- 
larize the solution and radically compress the data in the 
large-scale structure regime inS IV , and in the individual 
dark matter halo regime in § |y|. We discuss these results 
in§0. 
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II. FORMALISM 

We begin by briefly reviewing general mapmaking tech- 
niques in § II A, establishing notation used throughout 
the paper. These techniques are then appl ied to the two- 
dimensional weak lensing shear in § II B for the recon- 
struction of convergence maps. The latter follows the 
well-known Kaiser & Squires algorithm but we pay 
special attention to the propagation of noise into the con- 
vergence reconstruction as that will play a central r ole in 
the three-dimensional mapping that follows. In § [I C , 
we discuss the generalization to multiple source planes. 



A. Mapmaking 

Mapmaking can be formulated in terms of the general 
inverse problem j^], where we seek an estimate s a of a 
signal vector s Q from a data vector that is a linear 
projection P& a of the signal plus measurement noise n^, 



P ba s Q + n 6 . 



(1) 



The projection matrix P;, Q has dimensions (nb,n a ), the 
number of elements in the data and signal vectors respec- 
tively and hence need not be square. Here and below 
subscripts are labels and not elements of the vectors and 
matrices; elements will be denoted by [Pf, a ]ij- We as- 
sume that both the signal and the noise have zero mean 
(s Q ) = (ri(,) = 0, that the signal and noise are uncorre- 



cted, 



0, and that the noise covariance is known, 



N 



bb 



(n fc n b ; 



The statistical properties of the signal, 



may or may not be known. 

The estimated signal s a is that which minimizes 



where 



X 2 + H, 



X 2 = (d 6 - Pftai^'N^di, - P &a s a ) 



(2) 



(3) 



(4) 



(5) 



and the penalty function H is a set of constraints and/or 
a regularization to choose among degenerate solutions. 

Minimizing \ 2 (with H = 0) returns the linear estima- 
tor 



where 



R„ 



(6) 



(7) 



and this simple reconstruction is well-posed as long as 
the product in square brackets is invertible. If Pb a itself 
is invertible then R„b = Pr 1 and the estimator becomes 



independent of both the signal and the noise. The errors 
in the reconstruction, 



Sa — Sa — [R-afcPfca ~ I]s Q + Raftl*/, . 



imply a new noise covariance 
N„.„. 



((^a s a ) (s a S Q 



(8) 



(9) 



which is independent of the signal. Note that minimizing 
X 2 is not the same as minimizing the reconstruction noise 
N aa since reconstruction errors are not penalized by \ 2 
in the noise-dominated regime. In fact, it minimizes N aa 
subject to the constraint R a /,Pfc a = I M. 

For a noisy reconstruction, prior knowledge of the sta- 
tistical properties of the signal can be used as to set a 
penalty function for a new estimator of the signal 



H 



xt s -l- 



(10) 



Minimization of Eqn. (||) then returns the Wiener filtered 
estimate of the signal s w = R^d;,, where 



R 



wb 



= fs 



aa 1 
t 

aa*- bi 



PLN^Pfcap'PLN^ 

N^p 1 , (11) 



which has the heuristic form of signal/ (signal+noise) and 
so suppresses the signal in the noise-dominated regime. 
The estimator s^, has the noise properties 



N„ 



[R,„/,P 



wb^ba 



-fR^N^R 
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b^bb^-wb ■ 



(12) 



The Wiener estimator may alternately be derived as that 
which minimizes N^, [^) . 

For a well-posed inverse problem the minimum-^ 2 and 
Wiener reconstructions are related by an invertible op- 
eration. We can view the result of the former as provid- 
ing a new data vector d a = s a having noise n Q with 
a covariance N aa . Then with the model of the data 
d a = I s^, + n a , Eqn. ( pi] ) for the Wiener filter returns 



P^ada — Pu;a^a 5 



where 



P"Uja S aa [S Q 



N„ 



(13) 



(14) 



Since this matrix is invertible, the two reconstructions 
are formally equivalent. 

Because the minimum-^ 2 method does not require 
prior knowledge, we choose it as the primary mapping 
technique if the inverse problem is well-posed. Alterna- 
tively, for an ill-posed inverse problem, where the matrix 
P^Nj^Pfca is not invertible, Wiener filtering can serve 
as the primary technique. 

As the Wiener example implies, secondary processing 
operations on the primary map can be viewed as simply 
another round of mapmaking. Any linear operation that 
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is invertible will retain the same information content as 
the original so long as the noise covariance is properly 
propagated. We will use this technique to go from the 
observed lensing shear to the convergence to the density 
field and finally to the Wiener, signal-to-noise, or point 
source filtered density fields. 



B. Lensing Observables 

The distortion of images due to weak gravitational 
lensing is described by the Jacobian matrix of the map- 
ping between the two-dimensional source and image 
planes (e.g. §) 



1 - k - 71 -72 

-72 1 - K + 71 



(15) 



where all components are functions of position on the sky 
n. 

The galaxy ellipticities form a noisy estimator of the 
shear components (71,72) which are in turn related to 
the convergence k by (e.g. J9|) 

[ 7l ±i 72 ](n) = -i /"dn'^ K (n'), (16) 

where 8 = |n — n'| denotes the length of the pixel sep- 
aration vector, and <p denotes its azimuthal angle in the 
coordinate system that defines the shear components. By 
constructing a map of k from the shear data, one com- 
presses the data set by a factor of two and also transforms 
the data into a form that is more conveniently related to 
the three-dimensional density field (see § III). 

Discretizing the sky into pixels returns Eqn. ([If]) in 
the form of the general mapmaking problem of Eqn. (|l]), 
where the signal s K has been linearly projected onto the 
data space with measurement noise added to form the 
data vector, d 7 . Explicitly, if one orders the data vector 
as 

d 7 = { 7 i(ni),7 2 (n 2 );... ; 7i(n„ pix ), 7 2 (n„ pix )} , (17) 



the model becomes 

d-y P 7K s K 7 n 7 , 
with the projection matrix 

_ Aj cos 20jj 



(18) 



T«J( 2 j-l)j 



Aj sin 2(p i 



(19) 



where the indices run over the pixels of area Aj and the 
angles are defined as averages over the pixel. 

If the noise in the shear data is dominated by the in- 
trinsic ellipticity of galaxies, it is given by 



7 2 

n • • = rn • . frms 



[n ga i]i 



(20) 



where n ga i is a vector containing the number of galaxies 
per pixel. Here 7 rms is the rms error from intrinsic ellip- 
ticities and measurement errors per galaxy. The intrinsic 
alignment of galaxies on small scales is a potential source 
of correlated noise 0. While we neglect its currently 
uncertain contribution here, the framework we establish 
can handle any source of noise provided its covariance is 
known. 

We have implicitly assumed here a single set of shear 
data per pixel. If the full data set includes multiple ob- 
servations of the shear of varying quality or even simply 
the unbinned individual galaxy estimates themselves, one 
merely extends the data vector and the mapmaking algo- 
rithm combines them with the appropriate noise weight- 
ing. Mapmaking can also test the validity of Eqn. ( |To| ) , or 
more properly the data model of Eqn. fll8|), through a re- 
construction of the complementary "B-mode" map. This 
procedure amounts to multiplying the kernel in Eqn. ( |T6| ) 
by i corresponding to a rotation of the shear vectors by 
45°. The reconstructed map should be consistent with 
noise. 

The estimator of k and its noise properties then follow 
from Eqns. (^|) and @ aside from two subtleties. The 
first subtlety involves the so-called mass-sheet degener- 
acy: adding a constant to the n signal in Eqn. ( |l6| ) yields 
no effect in the shear data. There is then a singular value 
associated with the inversion in Eqn. (0). In general, one 
way to handle unconstrained modes is to add them to the 
data vector and assign them zero value but a large noise 
variance |0|. The covariance matrix of Eqn. (|^) then 
properly accounts for the lack of information on these 
modes. For the mass-sheet degeneracy, one appends: a 
zero to the data vector in Eqn. (|l7|); a row to the projec- 
tion matrix with uniform elements 1 / n p j x ; and a diagonal 
entry to the noise matrix Eqn. ( pp| ) with a value substan- 
tially greater than the variance of k smoothed across the 
field size in any reasonable cosmology. 

This same procedure applies to the second subtlety. 
Since Eqn. ( |l6| ) represents a convolution, the discrete 
representation is formally ill-defined for non-contiguous 
regions of the data, including holes and edges of the finite 
field. The sharp fall-off of the convolution kernel implies 
that only neighboring regions will be affected. Again 
one can account for these problems by assigning to the 
unmeasured or contaminated regions zero signal but a 
substantially larger noise variance than either the signal 
or noise in the neighboring measured region. The map- 
ping procedure then propagates an appropriately large 
and correlated noise into the reconstruction. If the gaps 
in the data are comparable to the contiguous regions, 
then Wiener filtering should serve as the primary map- 
ping technique. 

In the limit of infinitesimal pixels and an infinite con- 
tiguous field both these subtleties disappear and 

P*„P 7 «-I. (21) 
So if the noise in the shear is also uncorrelated and sta- 
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tistically homogeneous, N 77 cx I and Eqn. (^) becomes 

s K -> P 7K d 7 , (22) 

which is the discrete form of the Kaiser & Squires 
result. Furthermore, the noise in the n reconstruction 
N KK — > N 77 . We will use this approximation for illustra- 
tion purposes. These limiting behaviors and their map- 
making implications are simplest to derive in the Fourier 
domain (e.g. Jl2|). 



three-dimensional density maps (e.g. &). We shall see 
that in practice true mapping requires a prohibitively 
high signal-to-noise ratio in the lensing observables for 
reasons fundamental to the lensing projection. We fo- 
cus here on the radial reconstruction in a single angular 
pixel since the full three-dimensional distribution may be 
constructed as a collection of such reconstructions. 



A. Radial Mapmaking 



C. Multiple Source Planes 

So far we have implicitly assumed a single source plane 
for the lensing observables. In reality the source galaxies 
will be broadly distributed around some median set by 
the depth of the survey. It is this fact that makes three- 
dimensional mapping possible. 

For definitencss we will typically assume a median red- 
shift z mo d = 1 and a functional form [13] 



dN dD , , . 
dz dz 



(23) 



where D is the comoving distance in the fiducial cosmo- 
logical model, and D* is set to reproduce the median red- 
shift. The normalization constant A is chosen to match 
the number density of faint galaxies on the sky. We will 
take n = 3.6 x 10 5 deg -2 and 7 rms = 0.3 for illustra- 
tion purposes; this represents an estimate of the usable 
galaxies and the shear noise per galaxy measured from a 
space-based platform (A. Refregier, private communica- 
tion) . 

To separate the source galaxies into redshift bins, 
galaxy redshifts with errors that are smaller than the bin 
size are needed. Without spectroscopic redshifts, the pre- 
cision will be limited by photometric techniques. We will 
take a minimum redshift bin of Az — 0.025 to test the 
potential of future surveys. Current surveys return pho- 
tometric redshifts with Az w 0.06 ||. We shall see that 
reconstruction noise due to the finite number of galax- 
ies per bin dominates before this resolution is reached 
so that the photometric redshift errors of even current 
surveys are unlikely to be the limiting source of error. 

The mapmaking technique for multiple source planes 
with independent noise, as is appropriate for the intrinsic 
cllipticity noise of Eqn. (|20|), is the trivial generalization 
of a single source plane. With correlated noise in the 
shear from systematic effects or intrinsic galaxy align- 
ments, one forms a data vector of all the observations 
and applies the same mapmaking algorithm to estimate 
k in source redshift bins with appropriately correlated 
noise. 



III. THREE-DIMENSIONAL MAPPING 

Given two-dimensional convergence maps in multiple 
source planes, it is in principle possible to reconstruct full 



We now take as the data vector d K in n z source redshift 
bins (in a given angular pixel) and assume that its noise 
properties N KK are defined by the reconstruction in § II B. 



The model for the convergence k is a radial projection of 
the three-dimensional density distribution (e.g. H) 



K ( z s) = 2 H 0^r, 



, dD (D s - 

dz — — 

dz D 



D)D 



(24) 



where D is the comoving distance in a flat universe, sub- 
script s denotes evaluation at the redshift of the source z s , 
and A = S/a is the density fluctuation with the growth 
rate in a matter-dominated universe scaled out. Dis- 
cretizing Eqn. in redshift bins returns the general 
equation of mapmaking (see Eqn. ([!])) 



PkASa + n„ 




D, 
D, 



i > D j 
-i < Dj 



(25) 



'(26) 



where SDj is the width of bin j, the distances are mea- 
sured to the center of the bins, and we have offset the 
source redshift bins by 1 so that the projection matrix 
is purely lower triangular. For notational simplicity we 
have assumed that the binning in the signal and data 
space is the same, but the generalization is straightfor- 
ward. 

The reconstructed density field is then given by the 
general mapmaking equation 



where 



Here 



R 



An 



N 



AA 



PIaN-. 1 



N 



AA 



vN-iP^]- 1 , 



(27) 



(28) 



(29) 



is the noise covariance of the estimator. 

As noted by Taylor fl, the reconstruction is in princi- 
ple well-posed and does not require regularization if the 
density field is to be recovered to the same redshift reso- 
lution and range as the convergence data. Our more gen- 
eral treatment accounts for inhomogeneities and correla- 
tion in the noise, and even gaps in the data (see § B B|). 
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More importantly, it returns the noise covariance of the 
estimator. We shall see in the next section that with- 
out knowledge of the noise covariance the reconstructed 
density field cannot be used for any practical purpose. 

The multipixel generalization of radial mapmaking 
concatenates the vectors for each pixel. For uncorre- 
cted noise in K-pixels, the result is simply the appli- 
cation of mapmaking pixel-by-pixel since the noise ma- 
trix is then block diagonal in the pixels. For correlated 
noise, the noise matrices in the multipixel generalization 
of Eqns. (p8[) and couple neighboring pixels. 



B. Noise Properties 

To get a feel for the properties of the reconstruction, 
consider an idealization with redshift bins that are equal 
in comoving width SD, and noise in k that is both uncor- 
rected and homogeneous. Then 



where 



and 



P kA = CM, 
N KK = a 2 I, 



C=-(H 6D) 2 Q m , 



Ma = 



'(«+i-i)&r 
o 



j < i + 1 
else . 



(30) 



(31) 



(32) 



With these simplifications the reconstruction matrix is 

1 3 = i-2,i , 



i 127-1 
[RA«]<j = - 4: 7 X { - 2 J = * - 1 , 



C 2i- 1 



(33) 



else . 



Note that the projection matrix is lower triangular, and 
the reconstruction matrix is tridiagonal. 

For i ^> 1 the reconstructed density field is essen- 
tially the finite difference approximation of the second 
derivative of the convergence data; this is the same sec- 
ond derivative seen in Taylor's continuous method 
To understand this result, consider the response in n to 
a density fluctuation in a single redshift bin (also see 
Fig. |b). As we move out in redshift, the k response is 
zero until we reach the density fluctuation. As we cross 
the fluctuation k undergoes a sudden "acceleration," and 
thereafter grows slowly. With perfect data, one would 
identify density fluctuations as regions where the second 
derivative of k is large. The problem is that the k re- 
sponse accelerates from zero and hence will be hidden by 
noise locally. This is the fundamental limitation of weak 
lensing tomography with real data. 

Taking finite differences of the data amplifies the noise 
and strongly correlates it between neighboring pixels. For 



i > 1 the noise matrix is 



[N A A]y = ^2 



j = i - 2, i + 2 , 
3 = i - 1, i + 1 , 
3 = i , 
else . 



(34) 



The noise covariance has a very particular form, which 
corresponds to a finite-difference approximation of a 
fourth derivative (see Fig. |l|b). That the noise is cor- 
related in this specific way will turn out to be crucial in 
extracting any information from the reconstruction. 

To see this, consider a toy model where the galaxies are 
equally distributed among n z redshift bins that extend to 

1 /2 

a cosmologically interesting distance Vtm H D max = 1. 
Then a 1 — n z ^ ms /n ga \ and Q^HqSD = n" 1 , and the 
rms noise per bin in the reconstruction is 



5/2 



"gal 



3.6 x 10 5 



-1/2 



0.3 



(35) 



where we have scaled the result with numbers from S II C 



for degree scale pixels. Even with these generous assump- 
tions, the signal-to-noise per bin in the reconstruction 
is generally small and only approaches unity for den- 
sity fluctuations that approach unity when averaged over 
these large volumes. 

Notice that the noise per bin scales as n z , which 
seems to suggest that increasing the radial resolution de- 
creases the total signal-to-noise ratio. But that would 
be true only if the noise were uncorrelated. Instead, the 
fact that the noise covariance has a very specific form 
that is not seen in realistic signals allows it to be filtered 
out. The blind reconstruction of primary mapmaking is 
therefore useful only as a first step in the process. To ex- 
tract information out of the map, one must regularize the 
reconstruction with a prior assumption about the signal 
that is to be recovered. 



IV. LARGE-SCALE STRUCTURE 

For scales greater than about 10' and the depths 
reached by modern surveys, the convergence field in a 
typical region of sky is dominated by large-scale struc- 
tures in the underlying dark matter density field 
Moreover, the fluctuations are in the linear to quasilin- 
ear regime where theoretical modeling can be expected 
to give a go od p rior assumption about the statistics of 



the field (§ fVA). In this limit, one can quantify and 



better represent the information contained in the noisy 
three-dimensional density map obtained from the pri- 
mary mapmaking of the previous section. We ap ply two 
well-known techniques: Wiener filtering (§ IVCj) to rep- 



re sent th e map itself, and the Karhunen-Loeve transform 
(§ IV D) whose eigenmodes encapsulate and expose the 
underlying information contained in the map. 
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A. Signal Matrix 

In the linear regime, the signal matrix Saa is related 
to the linear power spectrum as follows. The average 
density fluctuation within the ith redshift window Wi (x) 
becomes 



or 1 / d 3 a;W i (x) ( 5(x), 



(36) 



where <$(x) is the density fluctuation field. We assume 
that the windows are normalized so that J d 3 xWi = 1. 
The signal covariance of these density averages is 



[S 



AA 



cPx, / (PxjWiixiJWjixij) 



x((5(x l )<5(x J )) 



(37) 



Gi Gj 

a,i a,j 



d 3 k 
(2^)3 



W t {k)W*(k)P(k), 



where P(k) is the linear power spectrum today Gi — 
G(zi) is the linear growth rate of the density field and 
Wi(k) are the Fourier transforms of the windows. Note 
that in a matter-dominated universe G{z) = a and so A 
then has the interpretation of the density field extrapo- 
lated to the present in linear theory. 

For definiteness, we will take the windows to be a series 
of slices in redshift at comoving d istan ce Di and width 
SDi with a sky pixel radius S = \J A/it in radians in the 
small angle approximation and a flat spatial geometry: 



W*(k) = 2e lk « D 



t Smk ll SD l /2J 1 (k ± D i e s ) 
k\\SDi/2 k ± D l Q s 



(38) 



For pixels smaller than A ~ 1 deg 2 , the density field 
is in the mildly non-linear regime E| . Here the signal 
matrix should be calculated with a numerical simulation 
of structure but a simpler approximation here suffices to 
extract the rough scaling. We replace in Eqn. d38|) 



G.G.jPik) -> [Pn^mO^WMj)] 172 



(39) 



where -Pnl(^) obtained from P(k) through the scaling 
relations of . This approximation says that structures 
maintain the same coherence as in linear theory across 
redshift bins. 



-0.5 




FIG. 1: True radial density distribution (shaded) vs. (a) 
Wiener and (b) primary reconstructions, for a 4 deg 2 pixel 
and source redshift bins with Az = 0.025 extending to 
Zmax = 2.5. Thin lines correspond to the fiducial noise vari- 
ance (n = 3.6 x 10 s deg~ 2 and 7 rms = 0.3), and lines of 
increasing thickness correspond to noise variance a factor of 
10, 100, and 1000 times lower. 



is a Gaussian realization of the correlated signal vector 
since 



(s A s t A )=L s (gg t )L* =Saa 



(42) 



B. Gaussian Realization 

The signal matrix can then be used to make a Gaussian 
simulation of structure. Consider the Cholesky decom- 
position of the signal matrix 

Saa = L S I4 , (40) 

and a vector of independent Gaussian random numbers 
of unit variance g, i.e. (gg 4 ) = I. Then 

sa = L s g (41) 



Again, below the degree scale, the Gaussian approxima- 
tion begins to break down. However as tested in simula- 
tions it remains a reasonable approximation down to 10' 




We show a sample realization in Fig. [j]a (shaded) for a 
ACDM cosmology with with parameters fi c = 0.3, fib = 
0.05, n A = 0.65, h = 0.65, n = 1, 5 H = 4.2 x 10~ 5 
((Ts = 0.92), pixel area of 4 deg 2 and redshift binning of 
Az = 0.025. We will use this as the fiducial cosmology 
in the examples that follow. 
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C. Wiener Filter 

As discussed in § II A, Wiener filtering minimizes the 
reconstruction noise using prior knowledge of Saa> the 
covariance matrix of the signal. In Fig. [I], we compare 
the Wiener reconstruction with the primary map for a 
Gaussian realization of structure and several choices of 
the noise variance. The thinnest lines (smoothest for 
Wiener, noisiest for primary) correspond to the fiducial 
noise spec ifications (n — 3.6 x 10 5 deg -2 and 7 rms = 0.3; 
see § II C ) , and the noise variance decreases by factors of 
10 as the lines thicken. Note the hundredfold difference 
in noise scale between the Wiener and primary recon- 
structions. 

The primary reconstruction is far too noisy to recover 
a visual impression of the structure, even for wildly opti- 
mistic assumptions about the noise. However, the noise 
has a specific oscillatory structure arising from the noise 
covariance (see Eqn. (|34|)), which is neither completely 
random nor present in the true signal. The Wiener fil- 
ter uses the information in the noise covariance of the 
primary reconstruction to reveal the hidden signal. For 
the fiducial noise specification, the Wiener filtered map 
recovers the low order, long-wavelength features in the 
density field; it is not until a prohibitively low noise vari- 
ance is reached that fine-scale features of order the bin 
width are recovered. 

The Wiener reconstruction is useful in cases where a 
map with well-defined statistical properties is needed, for 
example for cross-correlation studies with luminous trac- 
ers of the dark matter. 



D. KL Transform 

In the low signal-to-noise regime, it is more quanti- 
tatively useful to express the data in terms of a new 
set of orthogonal basis functions that are rank ordered 
by their signal-to-noise ratio. Low signal-to-noise modes 
may be eliminated from the data set allowing a near loss- 
less compression of the data. The pixel representation of 
these modes then tells us their correspondence to the ra- 
dial density field. This is accomplished by the Karhuncn- 
Loeve transform, also known as the signal-to-noise eigen- 
mode technique (e.g. [pj|). 

Consider the generalized eigenmode problem, 



SAAe = eNAAe . 
With a Cholcsky decomposition 



N 



AA 



LnL 



(43) 



(44) 



the generalized eigenmode problem reduces to an ordi- 
nary one 



[L^SaaLn 4 ] [L^e] = e[L*,e] 



(45) 



The eigenvectors represent linear combination of the data 
cIa- If one composes the matrix with rows representing 




eigen 



FIG. 2: Signal-to-noise ratio in the Karhunen-Loeve eigen- 
modes with the sky-coverage scaled out, (e/Zsky) 1 ^ 2 , calcu- 
lated for 5' x 5' pixels with the fiducial noise variance and 
redshift binning. 



the eigenvectors 



R 



KA 



(46) 



the new representation of the data vector becomes 

sk = RkacIa • (47) 

The important property of the Karhunen-Loeve trans- 
form is that 

(srSr) = Rka[S A a + N A a]R K a 

= Skk + Nkk, (48) 

where the the covariance matrices satisfy the important 
condition 



Skk 
N K k 



rl. 
I. 



(49) 



such that the modes are uncorrelated separately in each. 
Furthermore, e quantifies the relative contributions of sig- 
nal and noise in the mode. 

In Fig. |^, we show the signal-to-noise ratio per eigen- 
mode e 1 / 2 . This ratio is scaled by the square root of the 
fraction / sky of the sky covered by a survey, to reflect the 
increase in the signal-to-noise for a statistical detection 
with independent pixels. Here we have assumed a pixel 
size of 5' x 5', which is sufficiently small to extract most 
of the information on large-scale structure. 

Note the steep decrease in the signal-to-noise ratio as 
a function of the eigenmode index: the first few eigen- 
modes contain most of the information. This fact allows 
a radical compression of the data, here from 100 bins to a 
handful. Nonetheless, even though the ratio is small for 
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FIG. 3: Signal-to-noise eigenmodes for the five largest eigen- 
values in Fig. ^. The eigenmodes are similar to low frequency 
Fourier modes in the radial direction. 



duces to the same issue. Errors in our cosmological as- 
sumptions in the projection matrix make the primary 
map not correspond precisely to a density reconstruc- 
tion. Fortunately the form of the projection matrix is 
similar in all cosmologies: a broad bell-shaped weighting 
that peaks halfway to the source distance. Again, the 
well-defined linear operations involved allow us to pre- 
dict the statistics of the primary map given a cosmology 
in spite of the fact that it does not strictly represent the 
density field. 

Of course, if the recovered cosmology differs greatly 
from the assumed one, then the signal-to-noise eigen- 
modes will become an inefficient representation of the 
data. The best solution to both problems is to iterate 
the analysis and converge on a fiducial cosmology that 
fits the data. 



V. INDIVIDUAL DARK MATTER HALOS 



essentially all of the higher eigenmodes on the scale of 
individual pixels (/sky ~ 2 x 1CP 7 ), a survey comprising 
4-400 deg 2 (/sky ~ 10 _4 -10~ 2 ) has more than enough 
signal for a statistical detection. 

To understand the information stored in the higher 
modes, we plot the first few eigenmodes in redshift in 
Fig. [j] renormalized to have unit norm. The first eigen- 
vector simply shows the overall lensing efficiency when 
integrated over the whole source distribution, i.e. it has 
a single peak at a distance halfway to the median redshift 
z = 1. The higher eigenvectors increase the number of 
nodes with the boundary conditions that the weight is 
negligible near the observer at z = and well beyond the 
median redshift. They are therefore the analogues of low 
order Fourier wavevectors in the radial direction. 

These low order modes are sensitive to the cosmologi- 
cal model itself in that their values depend on the growth 
rate of structure, the volume element of the pixel-redshift 
bins, and the distances in the lensing efficiency. Even a 
statistical measure of their rms amplitude can help con- 
strain cosmology, in particular the dark energy. 

A potential problem for this use of three-dimensional 
mapping is that a cosmology is assumed both in the 
Karhunen-Loeve decomposition and in the projection 
matrix itself Eqn. (p6|). The former problem is readily 
handled in that even if the a priori assumption of the 
signal matrix in Eqn. ( p8| ) is incorrect, the Karhunen- 
Loeve transform is a well-defined linear operation on the 
data. The modification comes about in the calculation 
of the covariance of the estimators in Eqn. (ft8|) . The 
covariance is still diagonal in the noise but need not be 
diagonal in the signal nor do its elements have the inter- 
pretation of signal-to-noise. Still, it is calculable for the 
purpose of model fitting and does not present a funda- 
mental problem. 

The second problem is apparently more subtle but re- 



Below I 0', the convergence field is dominated by indi- 
vidual structures, or dark matter halos, along the line-of- 
sight to the source galaxies. The abundance in redshift of 
such objects is well known to be exceedingly sensitive to 
the growth rate of structure and hence the dark energy 
in the universe ]l8| . Identification and mass measure- 
ment by lensing would be ideal because the association 
of luminous observables with the dark mass of the ha- 
los is always problematic. Indeed there may exist halos 
that are effectively dark pl| . However, the efficacy of 
a purely lensing-based study is severely compromised by 
projection effects poj , so three-dimensional mapping in 
principle holds the key to utilizing this fundamental test. 

In the discrete halo case, one also has a well-motivated 
prior to regularize the inversion. 



In 



VA 



we discuss a 



modified version of the maximum entropy method. This 
method is most useful in the intermediate regime where 
the signal contains individual objects embedded in the 
large-scale structure. We then describe point source reg- 
ularization, which is the best method when there is good 
reason to believe that effectively all of the structures are 
well localized and the large-scale structure component 
can be ignored. 



A. Maximum Entropy Method 

The maximum entropy method (MEM) is widely ap- 
plied in situations where a noisy image is assumed to 
contain both discrete objects and a diffuse component 
(e.g. [ pT| ) and has been applied to two-dimensional weak 
lensing data [^2| . In this case the noisy image is the pri- 
mary map dA, the discrete objects are dark matter halos, 
and the diffuse component is the large-scale structure of 
the universe. MEM involves adding a penalty function to 
the mapmaking minimization of Eqn. (0) for the recovery 
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FIG. 4: MEM reconstruction with a halo of 5 x lO 14 ft -1 Ar 
added to a 10' x 10' pixel realization of large-scale structure 
(shaded) with the fiducial noise variance and redshift bin- 
ning (see Fig. Q). High (thick line, <5 t hr = 10) and low (thin 
line, <5 t h r = 3) threshold cases are shown. The MEM solution 
makes a discrete transition between these stable solutions. 



of an MEM filtered signal se, 



H = A£l[s E ] i ln/[ SE ] i 



(50) 



where the intensity / is some functional of se that we 
require to be positive. The Lagrange multiplier A trades 
off between minimizing x 2 and regularizing the solution. 
It is here chosen to give x 2 = v i the degrees of freedom. 
The main feature of MEM is that while it prefers a uni- 
form solution Ii =const., it does not additionally penal- 
ize bin-to-bin fluctuations as polynomial regularization 
would. Hence it allows solutions with discrete objects 
that occupy only one bin. 

To apply MEM to our mass reconstruction we need to 
choose I. A natural choice would be 1 + aj[sE]j since 
the density cannot fluctuate to negative values. However 
this prescription would still strongly disfavor placing all 
of the mass in a single redshift bin. To allow us more 
freedom in the regularization let us take a more general 
form 



I[x] = 1 + x, 
x = (5 t hr arctan 



a-i [se]. 



(51) 



The parameter 5thx places a threshold density above 
which MEM no longer penalizes the reconstruction; it 
returns the natural choice when (5 t hr - * co- 
in Fig. [| we show an example of the technique for the 
case of a halo embedded in large scale structure. For a 
high <5thr compared with the true density contrast of the 
halo, MEM seeks to regularize the reconstructed density 
in pixels and returns a smooth solution. For a <5thr that 



is lower, MEM places a density spike at the right posi- 
tion but the wrong amplitude. It favors a solution where 
the neighboring bins are all underdense since there is no 
penalty for further adding to the height of the spike. 

A fundamental drawback of MEM is the difficulty in 
assessing the errors in the reconstruction, i.e. the reality 
of the discrete objects MEM finds. Note that by con- 
struction both solutions in Fig. || have exactly the same 
X 2 , and the radical change in character of the solution 
is driven by the prior assumption of Sthr which sets the 
likelihood of having a comparable density spike in the 
solution. Still, MEM can identify interesting regions in 
the data for further study, perhaps with the point source 
method below. Conversely, it provides a useful cross 
check on the robustness of the point source solutions be- 
low. 



B. Point Source Method 

In regions that are known to be atypical of large-scale 
structure — either as flagged by the MEM reconstruc- 
tion or simply because the signal in k is much too large 
to be generated by large-scale structure — it is reason- 
able to assume as a prior that the density field is domi- 
nated by a collection of discrete massive objects. Ironi- 
cally, this form of anti-regularization of the density field 
in redshift bins is itself the most extreme regularization 
of the ones considered here, i.e. it has the least number 
of allowed degrees of freedom. The single-object form of 
this technique has been applied to data by || and yields 
impressively precise predictions of the redshift or radial 
location. Whether the predictions are accurate, however, 
depends on the validity of the single-object assumption 
and on the regularization criteria more generally. 

Let us state the criteria in a more general form. De- 
fine the penalty function on a point-source-regularized 
reconstruction sh as 



H = 



no < Tin 
oo no > n-u 



(52) 



where no is the number of redshift bins occupied by 
a density fluctuation, and nn is a prior assumption of 
the number of discrete objects (halos) the reconstruction 
should have. In other words, the minimization of Eqn. (Q) 
is strictly over the position and density amplitude of tir 
objects. The danger in this method of course is that it 
will return a best fit for the njj objects even if the solu- 
tion is in fact a smooth distribution or composed of some 
other number of objects. 

A well-defined procedure that makes minimal use of 
prior information is to identify sky pixels like to contain 
one or more massive objects (as described above), and to 
perform a sequence of minimizations with njj = 1,2,..., 
stopping when the x 2 does not improve significantly. 

In Fig. |^a, we show an example with two very massive 
cluster-sized halos at z = 0.19 and z = 0.34. (The masses 
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FIG. 5: Point source method, (a) Density reconstruction for 
two halos of mass 10 15 /i -1 Mq each added to a 5' x 5' pixel 
realization of large-scale structure (shaded), with the fiducial 
noise variance and binning (see Fig. Q). Priors of nn = 1 
(thin) and nn = 2 halos (thick) are compared, (b) Original 
data in the convergence « compared with the reconstructions. 



are enclosed entirely within the 5' x 5' pixels and their re- 
spective redshift bins.) The fit assuming one point source 
returns \ 2 — 96 for 100 redshift bins and two parameters 
— a perfectly good fit. The fit yields a redshift constraint 
z = 0.239 ± 0.006 that is remarkably precise, but wrong. 
Going to njj = 2 does recover two objects in the proper 
locations, with \ 2 — 91 for two fewer degrees of freedom. 
In Fig. ||b we show the implied k fields plotted against 
the original data. The residuals for the one-object fit 
show coherent structure near the true halo locations and 
so the improvement in x 2 is significant. Still, this exam- 
ple warns against blindly interpreting the formal errors of 
the fit. It is actually an optimistic example because the 
large masses and redshift separation (Az = 0.15) yield a 
signal much larger and much better separated than ex- 
pected for real weak lensing measurements. 

To better distinguish between close alternatives one 
could fold in prior information. For example, one could 
use the theoretically well-understood abundance of mas- 
sive halos to determine the relative likelihood of the solu- 



tions given the recovered masses of the objects, e.g. two 
halos of 5 x 10 15 h^ 1 Mq may be favored over one halo 
of 10 16 Mq due to the predicted exponential sup- 
pression in the number density of high mass halos. One 
could also use model profiles to create matched filters 
across smaller pixels that resolve the halo. Finally, prior 
information from photometric redshifts of galaxies likely 
to be members of the cluster (s) could decide between 
competing solutions j5j. 

VI. DISCUSSION 

We have shown that the evolution of the shear field 
in source redshift contains large-scale information about 
the distribution of dark matter, statistical information 
about its fluctuations on smaller scales, and redshift lo- 
calization information for massive dark matter halos on 
the smallest scales where the signal is large. This infor- 
mation is hidden in the noise of a direct reconstruction, 
and its extraction requires mild prior assumptions about 
the statistical properties of the density field. We have ar- 
gued for an approach that begins with a lossless direct or 
primary reconstruction that is followed by regularization 
by a prior that is appropriate for the information that is 
to be extracted. 

In the large to intermediate scale regime, the informa- 
tion content can be distilled into signal-to-noise or KL 
eigenmodes which efficiently compress the data by a fac- 
tor of 10 or more. These low-order modes probe the 
slow evolution of the statistics of the density field and 
are well suited to studying the properties of the dark 
energy. Tomographic sensitivity to the dark energy has 
been previously noted in the two-point correlation of the 
shear through the improvement of projected measures of 
the dark energy density and equation of state for future 
surveys p3[ . The KL eigenmode decomposition retains 
information from the higher order correlations in the field 
and also establishes a more direct, non-parametric 
quantification of the information contained in the data. 
A full study of the cosmological implications is beyond 
the scope of this paper, but we believe that it will be a 
promising approach for the future. 

Wiener filtering in the large-scale regime returns large 
scale maps of the density field with well-defined statisti- 
cal properties. These should be useful in cross-correlation 
studies with luminous, biased tracers of the dark matter 
such as galaxies and galaxy clusters |2f| . With informa- 
tion on the radial dimension, information on the evolu- 
tion of the bias can be recovered which in turn constrains 
the tracers' formation and evolution. 

In the individual halo regime, tomographic techniques 
have already been successfully applied to data || . With 
a well-motivated prior on the number of discrete halos 
along the line of sight, the reconstruction can yield ex- 
cellent localization of the object (s), in principle to a pre- 
cision that is better than that in the source redshifts 
themselves. However the accuracy is compromised by 
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an incorrect assumption of the number of objects. Here 
we advocate a combined approach of adding discrete ob- 
jects to the fit, regularizing by maximum entropy, and 
employing prior information and followup. 

While it is unfortunate that these prior assumptions 
are necessary for extracting information from three- 
dimensional reconstructions of the density field, they 
are generally well-motivated and testable. Gravitational 
lensing therefore remains our most direct, assumption- 
free means of probing the distribution of the dark mat- 



ter. 
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